Fisheries induced shifts in sea snake assemblages along the Konkan coast of India

Supplementary materials B: Analysis and Figures

options(max.print="75")
knitr::opts_chunk$set(echo = T, error = T, warning = F, message = F, tidy = T, cache = T,
                      dpi = 300, fig.align = "center", eval = T)

#Required libraries
suppressPackageStartupMessages(library(plyr))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(RColorBrewer))

#Importing sea snake bycatch data

snakes = read.csv("./Data/Rao et al_data_sea snake bycatch.csv")


# Fixing tow duration for 2018 - 2019

snakes <- snakes%>%
  mutate(Month = factor(months(dmy(Date), abbreviate = T), levels = month.abb))%>%
  group_by(Year, Gear.Type)%>%
  mutate(n.hauls = ifelse(is.na(No..of.Hauls), #filling in missing data
                         mean(No..of.Hauls, na.rm = T), No..of.Hauls),
         haul.time = ifelse(is.na(Average.Haul.Duration..Hours.), 
                         mean(Average.Haul.Duration..Hours., na.rm = T),
                         Average.Haul.Duration..Hours.))%>%
  ungroup()%>%
  mutate(Tow.duration.hours = 
           ifelse(is.na(Tow.duration.hours), # only for TD not recorded
                    ifelse(Gear.Type == "GillNet", #for gillnets
                            (haul.time+1.75), #soak time correction
                           ifelse(Gear.Type == "Trawler", #for trawlers
                                  n.hauls*haul.time, # total tow duration
                                  Tow.duration.hours) #return unchanged
                           ), #return unchanged
                    Tow.duration.hours)# return same for TD recorded
         )

#setting ggplot theme

theme_set(theme_bw()+
            theme(legend.position = "top",
                  text = element_text(size = 15)
                  )
          )

# function calculating mcfadden's pseudo r2 for glms

mcfadden <- function(x){
  r2 <- 1 - (x$deviance/x$null.deviance)
}

# function for calculating effect size for Z - test

proptest.r <- function(x, y){
  r <- x$statistic/sum(y$N)
}

Composition of sea snakes in bycatch

Sampling effort

We calculated the number of trips sampled and fishing effort per trip (Tow duration in hours) to determine sampling effort per month of each year of sampling.

#Number of boats sampled accross years

snakes%>%
  group_by(Year, Month)%>%
  distinct(Date, Boat.Name, .keep_all = T)%>%
  #filter(!Boat.Name == "")%>%
  count(name = "n.trips")%>%
  spread(Month, n.trips)
Year Jan Feb Mar Apr May Oct Nov Dec
2016 2 8 8 3 6 NA 4 10
2017 12 19 15 NA NA NA 1 6
2018 91 56 59 15 22 18 41 44
# Fishing effort accross years

snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  group_by(Gear.Type, Year, Month)%>%
  distinct(Date, Boat.Name, .keep_all = T)%>%
  summarise(sampling.effort = round(sum(Tow.duration.hours, na.rm = T), 2))%>%
  spread(Month, sampling.effort, fill = c(0))
Gear.Type Year Jan Feb Mar Apr May Oct Nov Dec
GillNet 2016 9.00 29.00 24.00 7.00 21.00 0.00 10.00 36.50
GillNet 2017 3.00 27.00 9.00 0.00 0.00 0.00 0.00 0.00
GillNet 2018 206.91 112.87 128.63 25.88 51.63 43.51 62.26 109.38
Trawler 2016 0.00 8.00 9.00 9.00 8.00 0.00 0.00 0.00
Trawler 2017 84.00 105.00 52.50 0.00 0.00 0.00 0.00 15.00
Trawler 2018 49.00 50.91 30.50 11.75 0.00 0.00 11.00 5.50
snakes%>%
  filter(Gear.Type == "Trawler" | Gear.Type == "GillNet")%>%
  group_by(Gear.Type)%>%
  distinct(Date, Boat.Name, .keep_all = T)%>%
  summarise(n.trips = n(),
            sampling.effort = round(sum(Tow.duration.hours, na.rm = T), 2),
            n.boats = length(unique(Boat.Name)))
Gear.Type n.trips sampling.effort n.boats
GillNet 329 916.57 136
Trawler 65 449.16 6

Sampling was not conducted trip wise for 2016 and most of 2017. Hence, effort data is unavailable.

Assemblage composition

We computed the abundance of each species of snake encountered over the whole sampling duration and across years

snakes%>%
  group_by(Species, Year)%>%
  count(name = "Abundance")%>%
  spread(Year, Abundance)
Species 2016 2017 2018
Acrochordus granulatus NA 1 NA
Hydrophis curtus 75 38 123
Hydrophis cyanocinctus NA NA 3
Hydrophis schistosus 178 215 521
Hydrophis viperinus NA 1 1

Varition in catch per unity effort (CPUE) of sea snakes across years in each gear

We considered only the two dominant species for in depth analysis. Catch Per Unit Effort (CPUE) for each species in each sampled trip was calculated as ‘Abundance of snakes encountered over Total Tow Duration or Soak Time’. We tested the effect of Tow duration on CPUE of each species using a log – linear regression. We compared the catch rates (CPUE) of each species across and within gears using log – linear regression.

CPUE.grsp <- snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus",
         Boat.Name != "")%>% # keeping dominant species
  droplevels()%>%
  group_by(Year, Gear.Type, Date, Boat.Name, Species)%>%
  summarise(Abn = n(), # abundance of sea snakes per trip
            fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort, mean for 2016 - 17
  complete(nesting(Year, Gear.Type, Date, Boat.Name, fe), 
           Species, fill = list(Abn = 0))%>% # completing species x trip combinations
  mutate(cpue = Abn/fe)%>%
  group_by(Boat.Name)%>%
  mutate(n = n())%>%
  ungroup()

CPUE.grsp%>%
  group_by(Gear.Type, Species)%>%
  summarise(Mean.CPUE = mean(cpue, na.rm = T))%>%
  spread(Species, Mean.CPUE)
Gear.Type Hydrophis curtus Hydrophis schistosus
GillNet 0.4243868 0.4647042
Trawler 1.1981818 1.7018479
CPUE.grsp%>%
  group_by(Gear.Type,Species)%>%
  nest()%>%
  mutate(shapiro = map(data, ~shapiro.test(.$cpue)),
         sumry = map(shapiro, broom::tidy))%>%
  select(sumry)%>%
  unnest()
Gear.Type Species statistic p.value method
GillNet Hydrophis schistosus 0.4785970 0.0000000 Shapiro-Wilk normality test
GillNet Hydrophis curtus 0.7492041 0.0000001 Shapiro-Wilk normality test
Trawler Hydrophis curtus 0.6313492 0.0001311 Shapiro-Wilk normality test
Trawler Hydrophis schistosus 0.7890101 0.0003384 Shapiro-Wilk normality test
# testing difference in CPUE between trawlers and gillnets

library(lme4)
library(broom.mixed)

CPUE.grsp%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~lmer(cpue ~ Gear.Type + (1|Boat.Name), data = .)),
         sumry = map(mod, broom.mixed::tidy))%>%
  select(sumry)%>%
  unnest()
Species effect group term estimate std.error statistic
Hydrophis schistosus fixed NA (Intercept) 0.4640419 0.0435137 10.664273
Hydrophis schistosus fixed NA Gear.TypeTrawler 1.2257846 0.1431319 8.564023
Hydrophis schistosus ran_pars Boat.Name sd__(Intercept) 0.0555275 NA NA
Hydrophis schistosus ran_pars Residual sd__Observation 0.6191548 NA NA
Hydrophis curtus fixed NA (Intercept) 0.4206151 0.1105983 3.803089
Hydrophis curtus fixed NA Gear.TypeTrawler 1.0609298 0.3247348 3.267065
Hydrophis curtus ran_pars Boat.Name sd__(Intercept) 0.4440595 NA NA
Hydrophis curtus ran_pars Residual sd__Observation 0.5613696 NA NA
CPUE.grsp%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~lmer(cpue ~ Gear.Type + (1|Boat.Name), data = .)),
         anva = map(mod, car::Anova),
         sumry = map(anva, broom::tidy))%>%
  select(sumry)%>%
  unnest()
Species term statistic df p.value
Hydrophis schistosus Gear.Type 73.34249 1 0.0000000
Hydrophis curtus Gear.Type 10.67372 1 0.0010867
# testing difference in CPUE between Species

CPUE.grsp%>%
  group_by(Gear.Type)%>%
  nest()%>%
  mutate(mod = map(data, ~lmer(cpue ~ Species + (1|Boat.Name), data = .)),
         sumry = map(mod, broom::tidy))%>%
  select(sumry)%>%
  unnest()
Gear.Type effect group term estimate std.error statistic
GillNet fixed NA (Intercept) 0.4244221 0.0437981 9.6904140
GillNet fixed NA SpeciesHydrophis schistosus 0.0402291 0.0485110 0.8292770
GillNet ran_pars Boat.Name sd__(Intercept) 0.0082453 NA NA
GillNet ran_pars Residual sd__Observation 0.3000941 NA NA
Trawler fixed NA (Intercept) 1.6265715 0.7494314 2.1704076
Trawler fixed NA SpeciesHydrophis schistosus -0.3410434 0.6845323 -0.4982138
Trawler ran_pars Boat.Name sd__(Intercept) 1.1674491 NA NA
Trawler ran_pars Residual sd__Observation 1.5111595 NA NA
CPUE.grsp%>%
  group_by(Gear.Type)%>%
  nest()%>%
  mutate(mod = map(data, ~lmer(cpue ~ Species + (1|Boat.Name), data = .)),
         anva = map(mod, car::Anova),
         sumry = map(anva, broom::tidy))%>%
  select(sumry)%>%
  unnest()
Gear.Type term statistic df p.value
GillNet Species 0.6877004 1 0.4069477
Trawler Species 0.2482170 1 0.6183334
#Plotting 

snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  droplevels()%>%
  group_by(Year, Gear.Type, Date, Boat.Name, Species)%>%
  summarise(Abn = n(), # abundance of sea snakes per trip
            fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort, mean for 2016 - 17
  complete(nesting(Year, Gear.Type, Date, Boat.Name, fe), 
           Species, fill = list(Abn = 0))%>% # completing species x trip combinations
  group_by(Year, Gear.Type, Species)%>%
  summarise(CPUE = mean(Abn/fe, na.rm = T), # Mean CPUE per trip
            sd = sd(Abn/fe, na.rm = T))%>% 
  ggplot(aes(Year, CPUE, shape = Species, group = Species))+
  geom_point(size = 3)+
  geom_line(aes(linetype = Species), size = 1)+
  scale_x_continuous(breaks = seq(2016, 2018, 1))+
  labs(y = "Mean CPUE")+
  facet_wrap(~Gear.Type)+
  theme(legend.text = element_text(face = "italic"))

Seasonal variation in CPUE

We tested for seasonality in bycatch trends for each species in each gear using log – linear regressions with months as predictor variables.

# testing seasonality in CPUE by gear

CPUE.month <- snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  droplevels()%>%
  group_by(Year, Month, Gear.Type, Date, Boat.Name, Species)%>%
  summarise(Abn = n(), # abundance of sea snakes per trip
            fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort, mean for 2016 - 17
  complete(nesting(Year, Month, Gear.Type, Date, Boat.Name, fe),  
           Species, fill = list(Abn = 0))%>% # completing species x trip combinations
  mutate(CPUE = Abn/fe) # Mean CPUE per trip

            
CPUE.month%>%
  ungroup()%>%
  group_by(Gear.Type, Species)%>%
  nest()%>%
  mutate(mod = map(data, ~lmer(CPUE ~ Month + (1|Boat.Name), data = .)),
         sumry = map(mod, broom::tidy),
         tst = map(mod, car::Anova))%>%
  select(sumry)%>%
  unnest()%>%
  select(term, estimate)%>%
  spread(term, estimate, drop = T)
Gear.Type Species (Intercept) MonthApr MonthDec MonthFeb MonthMar MonthMay MonthNov MonthOct sd__(Intercept) sd__Observation
GillNet Hydrophis curtus 0.3474585 1.2261151 0.1051078 0.0721653 1.6615975 NA 0.0344719 0.0343597 0.0000000 0.8967478
GillNet Hydrophis schistosus 0.4313831 0.0670926 0.0095525 0.5993502 0.0363549 0.0028650 -0.0718391 -0.0360249 0.0949908 0.5935458
Trawler Hydrophis curtus 0.3366162 -0.1013220 0.4050505 -0.0995791 -0.1143939 NA 3.8300505 NA 0.0000000 0.6425533
Trawler Hydrophis schistosus 0.7202427 1.6230328 0.9464240 0.3707643 0.4578930 0.1702764 2.0797573 NA 0.4216842 1.1841602
CPUE.month%>%
  ungroup()%>%
  group_by(Gear.Type, Species)%>%
  nest()%>%
  mutate(mod = map(data, ~lmer(CPUE ~ Month + (1|Boat.Name), data = .)),
         sumry = map(mod, broom::tidy),
         tst = map(mod, car::Anova))%>%
  select(tst)%>%
  unnest()
Gear.Type Species Chisq Df Pr(>Chisq)
GillNet Hydrophis schistosus 30.321943 7 0.0000829
Trawler Hydrophis schistosus 8.054082 6 0.2341679
GillNet Hydrophis curtus 23.811045 6 0.0005657
Trawler Hydrophis curtus 67.330098 5 0.0000000
#plotting

snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  droplevels()%>%
  mutate(Month = factor(Month, levels = month.abb))%>%
  group_by(Year, Month, Gear.Type, Date, Boat.Name, Species)%>%
  summarise(Abn = n(), # abundance of sea snakes per trip
            fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort, mean for 2016 - 17
  complete(nesting(Year, Gear.Type, Date, Boat.Name, fe), Month, 
           Species, fill = list(Abn = 0))%>% # completing species x trip combinations
  group_by(Month, Gear.Type, Species)%>%
  summarise(CPUE = mean(Abn/fe, na.rm = T), # Mean CPUE per trip
            sd = sd(Abn/fe, na.rm = T))%>%
  ggplot(aes(Month, CPUE, shape = Species,  group = Species))+
  geom_point(size = 3)+
  geom_line(aes(linetype = Species), size = 1)+
  scale_y_log10()+
  geom_vline(xintercept = "Jun", size = 1, linetype = "dotted")+
  geom_vline(xintercept = "Aug", size = 1, linetype = "dotted")+
  geom_text(aes(x = "Jul", y = 0.3, label = "Monsoon Ban"))+
  labs(y = "Mean CPUE (log)")+
  facet_wrap(~Gear.Type, ncol = 1)+
  theme(legend.text = element_text(face = "italic"))

Length distribution in bycatch

We determined selectivity of gears with respect to SVL of individuals using a t – test for each species in each gear.

#plotting

snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  ggplot(aes(Snout.to.Vent..cm., fill = Gear.Type))+
  geom_density(aes(linetype = Gear.Type),alpha = 0.3, size = 1)+
  labs(x = "Snout to vent length (cm)", y = "Kernel Density")+
  facet_wrap(~Species)+
  theme(strip.text = element_text(face = "italic"))

# summary SVL

snakes%>%
  filter(Species == "Hydrophis curtus" | Species == "Hydrophis schistosus")%>%
  group_by(Species)%>%
  skimr::skim(Snout.to.Vent..cm.)%>%
  skimr::yank("numeric")

Variable type: numeric

skim_variable Species n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Snout.to.Vent..cm. Hydrophis curtus 89 0.62 58.07 12.77 34.5 49.0 55.5 63.75 105.0 ▃▇▃▂▁
Snout.to.Vent..cm. Hydrophis schistosus 203 0.78 74.36 15.01 1.0 64.5 73.5 84.85 134.2 ▁▁▇▃▁
#Z - test for proportion

snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>%
  filter(!is.na(Snout.to.Vent..cm.))%>%
  droplevels()%>%
  select(Gear.Type, Species, Snout.to.Vent..cm.)%>%
  group_by(Species)%>%
  nest()%>%
  mutate(t.test = map(data, ~t.test(Snout.to.Vent..cm. ~ Gear.Type, data = .)),
         sumry = map(t.test, broom::tidy),
         d = map(data, ~lsr::cohensD(Snout.to.Vent..cm. ~ Gear.Type, data = .)))%>%
  unnest(sumry, d)%>%
  select(estimate:p.value, d)
Species estimate estimate1 estimate2 statistic p.value d
Hydrophis curtus -0.6316698 57.60896 58.24063 -0.2802908 0.7797036 0.0489142
Hydrophis schistosus -1.6728765 74.32615 75.99903 -1.4689341 0.1423528 0.1155450

Sex distibution in bycatch

We determined selectivity in terms of sex individuals using a two sample Z test for the proportion of females caught.

#Plotting

snakes%>%
  filter(Sex == "Male" | Sex == "Female")%>% # sex recorded
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  ggplot(aes(Species, fill = Sex))+
  geom_bar(position = "dodge", col = "black", width = 0.25)+
  scale_fill_brewer(palette = "Greys")+
  facet_wrap(~Gear.Type)+
  theme(axis.text.x = element_text(face = "italic"))

# Z - test for proportion

snakes%>%
  filter(Sex == "Male" | Sex == "Female")%>% # sex recorded
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>%
  group_by(Gear.Type, Species, Sex)%>%
  count()%>%
  spread(Sex, n)%>%
  mutate(N = Male+Female)%>%
  group_by(Species)%>%
  nest()%>%
  mutate(prop = map(data, ~prop.test(.$Female, .$N)),
         sumry = map(prop, broom::tidy),
         r = map2(prop, data, ~proptest.r(x = .x, y = .y)))%>%
  unnest(sumry, r)%>%
  select(estimate1, estimate2, p.value, conf.low, conf.high, r)
Species estimate1 estimate2 p.value conf.low conf.high r
Hydrophis curtus 0.5940594 0.2962963 0.0007665 0.1286188 0.4669074 0.0730372
Hydrophis schistosus 0.5082508 0.4444444 0.1848910 -0.0282980 0.1559108 0.0034468

Significantly more male H. curtus caught in trawler than gillnets

Effect of effort on composition of sea snakes

We tested the effect of Tow duration on CPUE of each species using a log – linear regression.

#Extracting CPUE data from bycatch database

CPUE = snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(!Boat.Name == "")%>% # removing points with no effort 
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  droplevels()%>%
  group_by(Gear.Type, Date, Boat.Name, Species)%>%
  summarise(Abn = n(), # abundance of sea snakes per trip
            fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort, mean for 2016 - 17
  complete(nesting(Gear.Type, Date, Boat.Name, fe), 
           Species, fill = list(Abn = 0))%>%
  mutate(CPUE = Abn/fe,
         log.Abn = log(Abn),
         log.CPUE = log(CPUE),# log CPUE because CPUE is 0 inflated
          log.fe = log(fe))# fe is 0 inflated

#Sample size

CPUE%>%
  group_by(Gear.Type)%>%
  distinct(Date, Boat.Name)%>%
  count()
Gear.Type n
GillNet 255
Trawler 27
#Plot

CPUE%>%
  ggplot(aes(log.fe, log.CPUE, shape = Species))+
  geom_jitter(size = 3)+
  geom_smooth(aes(linetype = Species), method = "lm", col = "black")+
  scale_color_brewer(palette = "Greys")+
  labs(y = "CPUE (log)", x = "Fishing Effort (log haul hours)")+
  facet_wrap(~ Gear.Type, scales = "free_x")+
  theme(legend.text = element_text(face = "italic"))

# testing 

CPUE%>%
  group_by(Gear.Type, Species)%>%
  nest()%>%
  mutate(mod = map(data, ~lmer(CPUE ~ fe + (1|Boat.Name), data = .)),
         sumry = map(mod, broom::tidy))%>%
  select(sumry)%>%
  unnest()
Gear.Type Species effect group term estimate std.error statistic
GillNet Hydrophis curtus fixed NA (Intercept) 0.7987359 0.1100043 7.260953
GillNet Hydrophis curtus fixed NA fe -0.1260171 0.0359139 -3.508871
GillNet Hydrophis curtus ran_pars Boat.Name sd__(Intercept) 0.0289787 NA NA
GillNet Hydrophis curtus ran_pars Residual sd__Observation 0.1792341 NA NA
GillNet Hydrophis schistosus fixed NA (Intercept) 0.6309144 0.1228337 5.136331
GillNet Hydrophis schistosus fixed NA fe -0.0596150 0.0433414 -1.375476
GillNet Hydrophis schistosus ran_pars Boat.Name sd__(Intercept) 0.0000000 NA NA
GillNet Hydrophis schistosus ran_pars Residual sd__Observation 0.3172430 NA NA
Trawler Hydrophis curtus fixed NA (Intercept) 3.4636588 1.7235253 2.009636
Trawler Hydrophis curtus fixed NA fe -0.3850708 0.2993008 -1.286568
Trawler Hydrophis curtus ran_pars Boat.Name sd__(Intercept) 1.3313764 NA NA
Trawler Hydrophis curtus ran_pars Residual sd__Observation 1.1400297 NA NA
Trawler Hydrophis schistosus fixed NA (Intercept) 3.9363791 0.9667292 4.071853
Trawler Hydrophis schistosus fixed NA fe -0.4062698 0.1380508 -2.942899
Trawler Hydrophis schistosus ran_pars Boat.Name sd__(Intercept) 0.7191601 NA NA
Trawler Hydrophis schistosus ran_pars Residual sd__Observation 1.4341989 NA NA
CPUE%>%
  group_by(Gear.Type, Species)%>%
  nest()%>%
  mutate(mod = map(data, ~lmer(CPUE ~ fe + (1|Boat.Name), data = .)),
         anva = map(mod, car::Anova),
         sumry = map(anva, broom::tidy))%>%
  select(sumry)%>%
  unnest()
Gear.Type Species term statistic df p.value
GillNet Hydrophis curtus fe 12.312177 1 0.0004500
GillNet Hydrophis schistosus fe 1.891933 1 0.1689840
Trawler Hydrophis curtus fe 1.655257 1 0.1982450
Trawler Hydrophis schistosus fe 8.660658 1 0.0032515

Mortality rates of sea snakes in bycatch by gear

We compared risk of mortalities of each species across and within gears using binomial regression with conditions at encounter as the response variable.

#Extracting mortality data from bycatch database

mortality <- snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  group_by(Gear.Type, Species)%>%
  summarise(n.alive = sum(Condition.at.encounter. == "Alive"),
            n.dead = sum(Condition.at.encounter. == "Dead"),
            mort.rate = (n.dead/(n.dead+n.alive))*100)

#table

mortality
Gear.Type Species n.alive n.dead mort.rate
GillNet Hydrophis curtus 97 34 25.954199
GillNet Hydrophis schistosus 448 34 7.053942
Trawler Hydrophis curtus 24 51 68.000000
Trawler Hydrophis schistosus 236 110 31.791908
#Plot

mortality%>%
  ggplot(aes(Species, mort.rate, fill = Gear.Type))+
  geom_col(position = "dodge", col = "black", width = 0.25)+
  scale_fill_brewer(palette = "Greys")+
  labs(x = "Species", y = "Mortality Rate (%)")+
  theme(axis.text.x = element_text(face = "italic"))

ggsave("./Figures/figure 3.tiff", device = "tiff", last_plot(), dpi = "print", 
       height = 4.5, width = 4.5)
# monthly mortalities

mortality.month <- snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  group_by(Month,Gear.Type, Species)%>%
  summarise(n.alive = sum(Condition.at.encounter. == "Alive"),
            n.dead = sum(Condition.at.encounter. == "Dead"),
            mort.rate = (n.dead/(n.dead+n.alive))*100)%>%
  ungroup()%>%
  mutate(Month = factor(Month, levels = month.abb))%>%
  complete(nesting(Gear.Type), Month, 
           Species, fill = list(mort.rate = 0)) # completing species x trip combinations
  
mortality.month%>%
  mutate(Month = factor(Month, levels = month.abb))%>%
  ggplot(aes(Month, mort.rate, shape = Species,  group = Species))+
  geom_point(size = 3)+
  geom_line(aes(linetype = Species), size = 1)+
  geom_vline(xintercept = "Jun", size = 1, linetype = "dotted")+
  geom_vline(xintercept = "Aug", size = 1, linetype = "dotted")+
  geom_text(aes(x = "Jul", y = 50, label = "Monsoon Ban"))+
  labs(y = "Mortality Rate (%)")+
  facet_wrap(~Gear.Type, ncol = 1)+
  theme(legend.text = element_text(face = "italic"))

ggsave("./Figures/figure 3b.tiff", device = "tiff", last_plot(), dpi = "print", 
       height = 4.5, width = 8)
#extracting mortality per trip from bycatch database

mort.grsp <- snakes%>%
  filter(Gear.Type == "GillNet" | Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  droplevels()%>%
  select(Gear.Type, Date, Boat.Name, Species, Condition.at.encounter.)%>%
  rename(Condition = Condition.at.encounter.)%>%
  mutate(Condition = as.factor(Condition))
# testing difference in mortality between trawlers and gillnets

mort.grsp%>%
  group_by(Species)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Condition ~ Gear.Type, data = ., family = binomial())),
         sumry = map(mod, broom::tidy),
         stat = map(mod, ~mcfadden(.)))%>%
  select(sumry, stat)%>%
  unnest()
Species term estimate std.error statistic p.value stat
Hydrophis schistosus (Intercept) -2.578433 0.1778870 -14.494775 0e+00 0.1131580
Hydrophis schistosus Gear.TypeTrawler 1.815081 0.2120660 8.559038 0e+00 0.1131580
Hydrophis curtus (Intercept) -1.048350 0.1993012 -5.260131 1e-07 0.1260767
Hydrophis curtus Gear.TypeTrawler 1.802122 0.3177978 5.670657 0e+00 0.1260767
# testing difference in CPUE between Species

mort.grsp%>%
  group_by(Gear.Type)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(Condition ~ Species, data = ., family = binomial())),
         sumry = map(mod, broom::tidy),
         stat = map(mod, ~mcfadden(.)))%>%
  select(sumry, stat)%>%
  unnest()
Gear.Type term estimate std.error statistic p.value stat
GillNet (Intercept) -1.0483505 0.1993014 -5.260126 0.0000001 0.0733598
GillNet SpeciesHydrophis schistosus -1.5300823 0.2671420 -5.727599 0.0000000 0.0733598
Trawler (Intercept) 0.7537718 0.2475368 3.045090 0.0023261 0.0596262
Trawler SpeciesHydrophis schistosus -1.5171232 0.2731349 -5.554484 0.0000000 0.0596262

Change in sea snake assemblages

Data on abundance of H. curtus and H. schistosus, overall sampling effort, number of trips was gathered from secondary published; Lobo (2004) and Padate et al (2009). A summary of the number of trips sampled and the sampling effort for each study in Table 2 of supplementary materials A. We used data from trawler bycatch alone to maintain consistency in our comparisons.

#importing data for past studies

past_studies <- read.csv("./Data/Rao et al_data_Past Studies.csv")

Sampling effort accross years

#summarising data from current study

effort.studies <- snakes%>%
  filter(Gear.Type == "Trawler")%>%
  group_by(Year)%>%
  distinct(Date, Boat.Name, .keep_all = T)%>%
  summarise(N.trips = n(),
            Total.fe = sum(Tow.duration.hours, na.rm = T))

#attching to previous studies

effort.studies <-  past_studies%>%
  select(Year, N.trips, Total.fe)%>%
  distinct(Year, .keep_all = T)%>%
  bind_rows(effort.studies)

#table

effort.studies
Year N.trips Total.fe
2002 194 385.00
2007 70 110.00
2016 5 34.00
2017 34 256.50
2018 26 158.66

Assemblage composition accross studies

We compared the relative abundance of H. curtus and mortality rates (proportion of individuals encountered dead) across studies using Z-Tests of proportion.

CPUE.studies <- snakes%>%
  filter(Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  droplevels()%>%
  group_by(Year, Date, Boat.Name, Species)%>%
  summarise(Abn = n(), # abundance of sea snakes per trip
            fe = mean(Tow.duration.hours, na.rm = T))%>% # fishing effort
  complete(nesting(Year, Date, Boat.Name, fe), 
           Species, fill = list(Abn = 0))%>% # completing species x trip combinations
  group_by(Year, Species)%>%
  summarise(Mean.CPUE = mean(Abn/fe, na.rm = T), # Mean CPUE per trip
            sd = sd(Abn/fe, na.rm = T),
            Total.Abn = sum(Abn, na.rm = T),
            Total.fe = sum(fe, na.rm = T))

CPUE.studies <- past_studies%>%
  select(Year, Species, Mean.CPUE, sd, Total.Abn, Total.fe)%>%
  bind_rows(CPUE.studies)
#Plotting

assemblage <- CPUE.studies%>%
  group_by(Year)%>%
  mutate(Total = sum(Total.Abn))%>%
  group_by(Year, Species)%>%
  summarise(rel.prop = Total.Abn/Total)%>%
  ggplot(aes(Year, rel.prop))+
  geom_point(aes(shape = Species), size = 3)+
  geom_line(aes(linetype = Species, group = Species), size = 1)+
  labs(y = "Relative Aundance in byactch")+
  theme(axis.text.x = element_blank())+
  theme(legend.text = element_text(face = "italic"),
        axis.title.x = element_blank())

# adding study periods to plot

study_periods <- read.csv("./Data/Rao et al_data_study_periods.csv")

# rectangle plot
period1 <-  ggplot(study_periods)+
  geom_rect(aes(xmin = start.year, xmax = end.year, ymax = 1, ymin = -1,
                fill = study),  col = "black", 
                 size = 0.1, height = 0.1)+
  geom_text(aes(x = median.year, y = 0, label = study), size = 3)+
  scale_fill_brewer(palette = "Greys", guide = F)+
  #scale_y_continuous(limits = c(0.5, 1.3))+
  labs(x = "Year")+
  theme_void()+
  theme(text = element_text(size = 15),
        axis.title.x = element_text(size = 15))

# arrow plot
period <-  ggplot(study_periods)+
  geom_errorbarh(aes(xmin = start.year, xmax = end.year, y = 1.75), 
                 size = 1)+
  geom_text(aes(x = median.year, y = 0.75, label = study), size = 3.5)+
  scale_y_continuous(limits = c(0, 2.25))+
  labs(x = "Year")+
  theme_void()+
  theme(text = element_text(size =15),
        axis.title.x = element_text(size = 15),
        axis.text.x = element_text())

# final plot
library(ggpubr)

ggpubr::ggarrange(assemblage, period, 
                  ncol = 1, 
                  heights = c(1, 0.2),
                  align = "v")

ggsave("./Figures/figure 4.tiff", device = "tiff", last_plot(), dpi = "print", 
       height = 4.5, width = 8)
# bar chart

CPUE.studies%>%
  mutate(Year = as.factor(Year))%>%
  group_by(Year)%>%
  mutate(Total = sum(Total.Abn))%>%
  group_by(Year, Species)%>%
  summarise(rel.prop = Total.Abn/Total)%>%
  ggplot(aes(Year, rel.prop))+
  geom_col(aes(fill = Species), col = "black", width = 0.5)+
  scale_fill_brewer(palette = "Greys")+
  labs(y = "Relative Aundance in byactch")+
  theme(legend.text = element_text(face = "italic"))

ggsave("./Figures/figure 4b.tiff", device = "tiff", last_plot(), dpi = "print", 
       height = 4.5, width = 8)
#Z - test for proportion

CPUE.studies%>%
  group_by(Year)%>%
  mutate(N = sum(Total.Abn))%>%
  filter(Species == "Hydrophis curtus")%>%
  ungroup()%>%
  nest()%>%
  mutate(proptest = map(data, ~prop.test(.$Total.Abn, .$N)),
         sumry = map(proptest, broom::tidy),
         r = map2(proptest, data, ~proptest.r(.x, .y)))%>%
  select(sumry,r)%>%
  unnest()%>%
  select(estimate1:p.value, r)
estimate1 estimate2 estimate3 estimate4 estimate5 statistic p.value r
0.8495575 0.2586207 0.1470588 0.1764706 0.1843318 281.6543 0 0.3995097

Difference in mortality of sea snakes in Trawlers across studies

#summarising data from current study

mortality.studies <- snakes%>%
  filter(Gear.Type == "Trawler")%>%
  filter(Species == "Hydrophis curtus" | 
           Species == "Hydrophis schistosus")%>% # keeping dominant species
  group_by(Year, Species)%>%
  summarise(n.alive = sum(Condition.at.encounter. == "Alive"),
            n.dead = sum(Condition.at.encounter. == "Dead"),
            mort.rate = (n.dead/(n.dead+n.alive)),
            Study = "Current")

# Attaching to previous study data

mortality.studies <- past_studies%>%
  select(Study, Year, Species, n.alive, n.dead, mort.rate)%>%
  bind_rows(mortality.studies)%>%
  mutate(N = n.alive + n.dead)%>%
  mutate(Study = factor(Study, levels = c("Lobo (2004)", "Padate et al. (2009)", "Current")))%>%
  group_by(Study, Species)%>%
  summarise(n.dead = sum(n.dead),
            N = sum(N),
            mort.rate = n.dead/N)
  
#table

mortality.studies
Study Species n.dead N mort.rate
Lobo (2004) Hydrophis curtus 132 192 0.6875000
Lobo (2004) Hydrophis schistosus 4 34 0.1176471
Padate et al. (2009) Hydrophis curtus 15 NA NA
Padate et al. (2009) Hydrophis schistosus 43 NA NA
Current Hydrophis curtus 51 75 0.6800000
Current Hydrophis schistosus 110 346 0.3179191
#plot

mortality.studies%>%
  ggplot(aes(Study, mort.rate, fill = Species))+
  geom_col(position = position_dodge(preserve = "single"), col = "black", width = 0.25)+
  scale_fill_brewer(palette = "Greys")+
  labs(x = "Year", y = "Mortality Rate (%)")+
  theme(legend.text = element_text(face = "italic"))

# test

mortality.studies%>%
  group_by(Species)%>%
  drop_na()%>%
  nest()%>%
  mutate(proptest = map(data, ~prop.test(.$n.dead, .$N)),
         #mod = map(data, ~lm(rel.prop ~ Year, data = .)),
         sumry = map(proptest, broom::tidy),
         r = map2(proptest, data, ~proptest.r(.x, .y)))%>%
  select(sumry,r)%>%
  unnest()%>%
  select(estimate1:p.value, r)
Species estimate1 estimate2 statistic p.value r
Hydrophis curtus 0.6875000 0.6800000 0.000000 1.0000000 0.0000000
Hydrophis schistosus 0.1176471 0.3179191 4.997571 0.0253829 0.0131515

Increase in fishing pressure in Konkan 2000 - 2020

In order to understand the change in fishing pressure in the region we gathered data on the number of vessels of each craft type (mechanised, motorised and non- motorised) for Maharashtra and Goa from CMFRI Marine Fisheries Census for 2005 and 2010; Fisheries Department for 2014, 2015, 2017 and 2018; Jena and George (2018) for 2016 and the Handbook on Fisheries Statistics for 2019. The CMFRI Marine Fisheries Census recorded all vessels across multiple harbours across each state. The other sources reported on registered vessels only. Change in the number of vessels of each craft type from 2005 to 2019 were tested using log – linear regression. We compared the number of vessels of each craft across Maharashtra and goa using a t – test.

#importing data on number of vessels

vessels <- read.csv("./Data/Rao et al_data_state-craft-gear.csv")
#plot

vessels%>%
  ggplot(aes(Year, n, shape = Craft))+
  geom_point(size = 3)+
  geom_line(aes(group = Craft, linetype = Craft), size = 1)+
  scale_y_sqrt()+
  labs(y = "No. of Vessels (sq. rt.)")+
  facet_wrap(~ State, ncol = 1)

ggsave("./Figures/figure 5.tiff", device = "tiff", last_plot(), dpi = "print", 
       height = 6, width = 8)
#testing change in number of vessels over years

vessels%>%
  group_by(State, Craft)%>%
  nest()%>%
  mutate(mod = map(data, ~glm(n ~ Year, data = ., family = poisson())),
         sumry = map(mod, broom::tidy),
         r2 = map(mod, ~mcfadden(.)))%>%
  select(sumry, r2)%>%
  unnest()%>%
  arrange(desc(r2))
State Craft term estimate std.error statistic p.value r2
Maharashtra Mechanised (Intercept) -69.9447237 1.3999812 -49.961190 0 0.7467118
Maharashtra Mechanised Year 0.0395470 0.0006949 56.913780 0 0.7467118
Goa Mechanised (Intercept) -115.4983459 4.2006942 -27.495062 0 0.5441047
Goa Mechanised Year 0.0611171 0.0020846 29.318823 0 0.5441047
Goa Non - motorised (Intercept) 86.2325473 8.6498363 9.969269 0 0.3073411
Goa Non - motorised Year -0.0399799 0.0042961 -9.306161 0 0.3073411
Goa Motorised (Intercept) 252.8613557 5.8641018 43.120219 0 0.2632877
Goa Motorised Year -0.1224781 0.0029154 -42.010095 0 0.2632877
Maharashtra Motorised (Intercept) 185.9700181 3.4680508 53.623788 0 0.1044864
Maharashtra Motorised Year -0.0886872 0.0017235 -51.458983 0 0.1044864
Maharashtra Non - motorised (Intercept) -33.6241328 2.2838299 -14.722696 0 0.0459403
Maharashtra Non - motorised Year 0.0209820 0.0011336 18.508824 0 0.0459403
#testing difference in number of vessels between Mah and Goa

vessels%>%
  group_by(Craft)%>%
  nest()%>%
  mutate(mod = map(data, ~t.test(n ~ State, data = .)),
         sumry = map(mod, broom::tidy),
         d = map(data, ~lsr::cohensD(n ~ State, data = .)))%>%
  select(sumry, d)%>%
  unnest()%>%
  arrange(desc(d))%>%
  select(estimate:p.value, d)
Craft estimate estimate1 estimate2 statistic p.value d
Mechanised -14568.29 2054.429 16622.714 -10.598781 0.0000233 5.6652869
Non - motorised -5365.00 304.625 5669.625 -6.387398 0.0003649 3.1936991
Motorised -1103.25 558.750 1662.000 -1.260113 0.2435439 0.6300565

Study site map

library(tmap)
library(osmdata)

study_sites <- read.csv("./Data/Rao et al_data_study_site.csv")

# geocoding study sites

library(ggmap)

study_sites <- study_sites%>%
  filter(Study != "2006-2008")%>%
  mutate(Name = as.character(Name))%>%
  mutate_geocode(Name)

study_sites
Study Name lon lat
2016 - 2018 Malvan 73.47105 16.06307
2002 - 2003 Chapora 73.74337 15.60212
2002 - 2003 Vasco da Gama, Mormugao 73.81133 15.39819
2002 - 2003 Malim Jetty 73.83445 15.50495
2002 - 2003 Betul, Ambelim 73.95637 15.14703
2016 - 2018 Vengurla 73.63890 15.85141
# converting to simplfe feature

library(sf)

study_sites <- st_as_sf(study_sites, coords = c("lon", "lat"))

st_crs(study_sites) <- 4326

# getting study extent

study_ext <- st_bbox(study_sites)

# getting osm map

library(osmdata)

land_raw <- opq(bbox = study_ext, timeout = 100)%>%
  add_osm_feature(key = 'admin_level', value = '4')%>%
  osmdata_sf()%>%
  unique_osmdata()
  
land <- land_raw$osm_multipolygons

st_crs(land) <- 4326

india <- read_sf("./Data/India-States.shp")

# plotting site map

ss_main <- tm_shape(land)+
  tm_polygons()+
  tm_text("name", remove.overlap = F)+
  tm_shape(india)+
  tm_polygons()+
tm_shape(study_sites, bbox = study_ext+c(-0.5,-.5,0.5,0.5), is.master = T)+
  tm_symbols(size = 1, shape = "Study", border.lwd = 2.5, border.col = "black", col = "gray")+
  tm_text("Name", size = 0.7,fontface = "bold",
          ymod = -0.7, xmod = -0.75, remove.overlap = T)+
  tm_graticules(n.x = 5)+
  tm_scale_bar(position = c("left", "top"), text.size = 0.5) + 
  tm_compass(position = c('right', 'top'), size = 5)+
  tm_style(style = "gray")+
  tm_layout(scale = 1.5, legend.position = c("left", "bottom"))

# making inset map

ext <- st_as_sfc(study_ext+c(-1.2,-1.2,1.2,1.2))

ss_inset <- tm_shape(india)+
  tm_polygons()+
tm_shape(ext)+
  tm_borders(lwd = 3)+
  tm_style("gray")+
  tm_layout(scale = 1.5)

# making final map

vp <- grid::viewport(x = 0.75, y = 0.15, width = 0.2, height = 0.3)

tmap_save(ss_main, insets_tm = ss_inset, insets_vp = vp, 
          "./Figures/Figure_1.tiff", height = 8, width = 8)